*
* $Id$
*
* $Log: stpair.F,v $
* Revision 1.1.1.1  2002/06/16 15:18:40  hristov
* Separate distribution  of Geant3
*
* Revision 1.1.1.1  1999/05/18 15:55:19  fca
* AliRoot sources
*
* Revision 1.1.1.1  1995/10/24 10:20:59  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.38  by  S.Giani
*-- Author :
      SUBROUTINE STPAIR
C
C *** STRANGE PARTICLE PAIR PRODUCTION ***
C *** NVE 14-MAR-1988 CERN GENEVA ***
C
C ORIGIN : H.FESEFELDT 16-DEC-1987
C
C THE SAME FORMULA FOR <K KB> VS AVAILABLE ENERGY
C                  AND <K Y>  VS AVAILABLE ENERGY
C FOR ALL REACTIONS.
C CHOOSE CHARGE COMBINATIONS K+ K- , K+ K0B, K0 K0B OR K0 K-
C                            K+ Y0, K0 Y+, K0 Y-
C FOR ANTIBARYON INDUCED REACTIONS HALF OF THE CROSS SECTIONS
C KB YB PAIRS ARE PRODUCED
C CHARGE IS NOT CONSERVED , NO EXPERIMENTAL DATA AVAILABLE FOR
C EXCLUSIVE REACTIONS, THEREFORE SOME AVERAGE BEHAVIOUR ASSUMED.
C THE RATIO L/SIGMA IS TAKEN AS 3:1 (FROM EXPERIMENTAL LOW ENERGY)
C
#include "geant321/s_defcom.inc"
C
      REAL KKB,KY
      DIMENSION KKB(9),KY(12),IPAKKB(2,9),IPAKY(2,12),IPAKYB(2,12)
      DIMENSION AVKKB(12),AVKY(12),AVNNB(12),AVRS(12)
      DIMENSION RNDM(1)
      DATA KKB/0.2500,0.3750,0.5000,0.5625,0.6250,0.6875,0.7500,
     *         0.8750,1.000/
      DATA KY /0.200,0.300,0.400,0.550,0.625,0.700,0.800,0.850,
     *         0.900,0.950,0.975,1.000/
      DATA IPAKKB/10,13,10,11,10,12,11,11,11,12,12,11,12,12,
     *            11,13,12,13/
      DATA IPAKY /18,10,18,11,18,12,20,10,20,11,20,12,21,10,
     *            21,11,21,12,22,10,22,11,22,12/
      DATA IPAKYB/19,13,19,12,19,11,23,13,23,12,23,11,24,13,
     *            24,12,24,11,25,13,25,12,25,11/
      DATA AVRS/3.,4.,5.,6.,7.,8.,9.,10.,20.,30.,40.,50./
      DATA AVKKB/0.0015,0.005,0.012,0.0285,0.0525,0.075,0.0975,
     *           0.123,0.28,0.398,0.495,0.573/
      DATA AVKY /0.005,0.03,0.064,0.095,0.115,0.13,0.145,0.155,
     *           0.20,0.205,0.210,0.212/
      DATA AVNNB/0.00001,0.0001,0.0006,0.0025,0.01,0.02,0.04,
     $           0.05,0.12,0.15,0.18,0.20/
C
      IF(IPA(3).LE.0) GO TO 9999
      IER(50)=IER(50)+1
      IPA1=ABS(IPA(1))
      IPA2=ABS(IPA(2))
C --- PROTECTION AGAINST ANNIHILATION PROCESSES ---
      IF ((IPA1 .EQ. 0) .OR. (IPA2 .EQ. 0)) GO TO 9999
      EAB=RS-ABS(RMASS(IPA1))-ABS(RMASS(IPA2))
      IF(EAB.LT.1.) GO TO 9999
C**
C** CHOOSE RANDOM REPLACEMENT OF PRODUCED KAONS (16.12.87)
      DO 111 I=1,60
      IF(IPA(I).EQ.0) GOTO 112
  111 CONTINUE
  112 I=I-3
      CALL GRNDM(RNDM,1)
      I3=3+IFIX(RNDM(1)*I)
  114 CALL GRNDM(RNDM,1)
      I4=3+IFIX(RNDM(1)*I)
      IF(I.EQ.1) I4=4
      IF(I3.EQ.I4) GOTO 114
C
C *** CHOOSE RANDOM REPLACEMENT OF PRODUCED KAONS (16.12.87) ***
C --- GET RS BIN ---
      DO 1 I=2,12
      IF (RS .LE. AVRS(I)) GO TO 2
 1    CONTINUE
      I1=11
      I2=12
      GO TO 3
C
 2    CONTINUE
      I1=I-1
      I2=I
C
C *** USE LINEAR INTERPOLATION OR EXTRAPOLATION BY Y=RC*X+B ***
 3    CONTINUE
      DXNVE=AVRS(I2)-AVRS(I1)
      DYNVE=LOG(AVKKB(I2))-LOG(AVKKB(I1))
      RCNVE=DYNVE/DXNVE
      BNVE=LOG(AVKKB(I1))-RCNVE*AVRS(I1)
      AVK=RCNVE*RS+BNVE
      DYNVE=LOG(AVKY(I2))-LOG(AVKY(I1))
      RCNVE=DYNVE/DXNVE
      BNVE=LOG(AVKY(I1))-RCNVE*AVRS(I1)
      AVY=RCNVE*RS+BNVE
      DYNVE=LOG(AVNNB(I2))-LOG(AVNNB(I1))
      RCNVE=DYNVE/DXNVE
      BNVE =LOG(AVNNB(I1))-RCNVE*AVRS(I1)
      AVN  =RCNVE*RS+BNVE
C
      AVK=EXP(AVK)
      AVY=EXP(AVY)
      AVN=EXP(AVN)
      IF(AVK+AVY+AVN.LE.0.) GOTO 9999
      IF(IPA1.LT.14) AVY=AVY/2.
      IF(IPA2.LT.14) AVY=0.
      AVY=AVY+AVK+AVN
      AVK=    AVK+AVN
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      IF(RAN.LT.AVN) GOTO 5
      IF(RAN.LT.AVK) GOTO 10
      IF(RAN.LT.AVY) GOTO 20
      GO TO 9999
    5 IF((EAB-2.).LT.0.) GO TO 9999
      CALL GRNDM(RNDM,1)
      IF(RNDM(1).LT.0.5) GO TO 6
      IPA(I3)=14
      IPA(I4)=15
      GOTO 30
    6 IPA(I3)=16
      IPA(I4)=17
      GOTO 30
   10 IF((EAB-1.).LT.0.) GO TO 9999
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      DO 11 I=1,9
      IF(RAN.LT.KKB(I)) GOTO 12
   11 CONTINUE
      GO TO 9999
   12 IPA(I3)=IPAKKB(1,I)
      IPA(I4)=IPAKKB(2,I)
      GOTO 30
   20 IF((EAB-1.6).LT.0.) GO TO 9999
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      DO 21 I=1,12
      IF(RAN.LT.KY(I)) GOTO 22
   21 CONTINUE
      GO TO 9999
   22 IF(IPA(1).LT.14) GOTO 23
      CALL GRNDM(RNDM,1)
      IF(RNDM(1).LT.0.5) GOTO 23
      IPA1=ABS(IPA(1))
      IPA(1)=IPAKY(1,I)
      IF(IPA1.EQ.15) GOTO 25
      IF(IPA1.EQ.17) GOTO 25
      IF(IPA1.EQ.19) GOTO 25
      IF(IPA1.GT.22) GOTO 25
      GOTO 24
   25 IPA(1)=IPAKYB(1,I)
      IPA(I3)=IPAKYB(2,I)
      GOTO 30
   23 IPA(2)=IPAKY(1,I)
   24 IPA(I3)=IPAKY(2,I)
C** CHECK THE AVAILABLE ENERGY
   30 EAB=RS
      IJ=0
      DO 31 I=1,60
      IF(IPA(I).EQ.0) GOTO 31
      IPA1=ABS(IPA(I))
      EAB=EAB-ABS(RMASS(IPA1))
      IJ=IJ+1
      IF(EAB.LT.0.) GOTO 35
   31 CONTINUE
      IF (NPRT(4)) WRITE(NEWBCD,1003) (IPA(J),J=1,IJ)
      GO TO 9999
   35 I=I-1
      L=I-1
      IF(L.LE.0) GO TO 9999
      DO 36 J=I,60
   36 IPA(J)=0
      IF (NPRT(4)) WRITE(NEWBCD,1002) (IPA(J),J=1,L)
C
 1002 FORMAT(' *STPAIR* KKB/KY PAIR PRODUCTION NOT ENOUGH ENERGY',
     $ ' REDUCE NUMBER OF PARTICLES ',2X,20I3)
 1003 FORMAT(' *STPAIR* KKB/KY PAIR PRODUCTION ENERGY SUFFICIENT',
     $ ' NUMBER OF PARTICLES ',2X,20I3)
C
 9999 CONTINUE
C
      END
